import numpy as np
import matplotlib.pyplot as plt

fig, (ax_diff, ax_aeff) = plt.subplots(1,2, figsize=(6,3))

# the differential limit and sensitivity
diff_file = np.genfromtxt("differential_limit_and_sensitivity.csv",
                          skip_header=1, delimiter=",",
                          names=["energy", "limit", "sensitivity"]
                          )
ax_diff.plot(diff_file["energy"], diff_file["sensitivity"], label="sensitivity")
ax_diff.plot(diff_file["energy"], diff_file["limit"], label="limit")
ax_diff.loglog()
ax_diff.legend()
ax_diff.set_xlabel("Neutrino Energy / GeV")
ax_diff.set_ylabel(r'E$^2 \cdot \Phi_{\mathrm{all-flavor}}^{\nu + \bar{\nu}}$ / (GeV cm$^{-2}$ sr$^{-1}$ s$^{-1}$)')


# the effective area
aeff_file = np.genfromtxt("effective_area.csv",
                          skip_header=1, delimiter=",",
                          names=["energy", "total", "nue", "numu", "nutau"]
                          )
ax_aeff.plot(aeff_file["energy"], aeff_file["total"], label="total")
ax_aeff.plot(aeff_file["energy"], aeff_file["nue"], '--', label=r"$\nu_e$")
ax_aeff.plot(aeff_file["energy"], aeff_file["numu"], '-.', label=r"$\nu_{\mu}$")
ax_aeff.plot(aeff_file["energy"], aeff_file["nutau"], ':', label=r"$\nu_{\tau}$")
ax_aeff.loglog()
ax_aeff.legend()
ax_aeff.set_xlabel("Neutrino Energy / GeV")
ax_aeff.set_ylabel(r'$A_{\mathrm{eff}} \, / \, \mathrm{m}^2$')

plt.tight_layout()
fig.savefig("demo.png")
